### Code for: Spatial-temporal clustering analysis of yaws on Lihir Island, Papua New Guinea
### to enhance planning and implementation of eradication programs
### Eric Q. Mooring, Oriol Mitja, and Megan B. Murray, 2018
### Code 1 of 3
### Spatial-temporal clustering of passively detected yaws cases at Lihir Medical Centre ###

# Ensure that the package "rsatscan" is installed.
# Ensure that SaTScan is installed.
# Make a directory called "SaTScan_Files" to put the input and output files for SaTScan

sslocdir<-"C:/Program Files/SaTScan" # Update name of directory where SaTScan is located

library(rsatscan)

### Prepare cases data ###

YawsCases<-read.csv(file = "YawsCases_LMC.csv")
YawsCases<-cbind(YawsCases,"NumCases" = rep(1,times=nrow(YawsCases))) # make row of 1s to be number of cases

total_episodes<-nrow(YawsCases)
YawsCasesConf<-YawsCases[YawsCases$Titer!="ND",] # Drop yaws disease episodes without a recorded RPR test for sensitivity analyses
total_episodes - nrow(YawsCasesConf) # drop 261

# SaTScan will automatically add the number of cases sharing a given attribute, so it's not necessary to tabulate them here

YawsCasesFile<-YawsCases[,c("Location","NumCases","Date")] # Make case file without covars; select and reorder columns for SatScan 
YawsCasesFileCovars<-YawsCases[,c("Location","NumCases","Date","AgeCat","Gender")] # Make case file with covars; select and reorder columns for SatScan 
YawsCasesFileCovarsConf<-YawsCasesConf[,c("Location","NumCases","Date","AgeCat","Gender")] # Make case file with covars, confirmed cases only; select and reorder columns for SatScan 

### Prepare "population" (all outpatient visits) data ###

# Load csv of all outpatient visits (visits from people originating from a non-Lihir village and duplicates have already been removed)
AllVisits<-read.csv(file = "./AllOutpatient_LMC.csv")

# Make pop table WITHOUT covars for Poisson analysis

AllVisitsTable<-table(AllVisits$Location, AllVisits$Date) # make table crosstabulating by date and village, only
AllVisitsFile<-as.data.frame(AllVisitsTable) # turn table into dataframe
colnames(AllVisitsFile)<-c("Location","Date","NumVisits") # relabel the columns of dataframe; SaTScan needs: location ID, time, population in this order

# Make pop table WITH covars for Poisson analysis

AllVisitsTableCovars<-table(AllVisits$Location, AllVisits$Date, AllVisits$AgeCat, AllVisits$Gender) # make table crosstabulating by date and village and age category and gender
AllVisitsFileCovars<-as.data.frame(AllVisitsTableCovars) # turn table into dataframe
AllVisitsFileCovars<-AllVisitsFileCovars[,c(1,2,5,3,4)] # reorder columns so number of visits ("population") is third column
colnames(AllVisitsFileCovars)<-c("Location","Date","NumVisits","AgeCat","Gender") # relabel the columns of dataframe; SaTScan needs: location ID, time, population, covars

### Prepare non-Euclidean neighbors file from list of distances between each pair of villages ###

VillageDistances<-read.csv(file="LihirVillageRoadDistances.csv") # load the village distances data
VillagesInPassive<-which((VillageDistances$FromVillage %in% levels(YawsCases$Location)) & (VillageDistances$ToVillage %in% levels(YawsCases$Location)))
VillageDistances<-VillageDistances[VillagesInPassive,]
NonEuclideanNeighbors<-matrix(data=VillageDistances$ToVillage,nrow=26,ncol=26,byrow=T) # make matrix, which has a column for each location
NonEuclideanNeighbors<-as.data.frame(NonEuclideanNeighbors) # save matrix as dataframe (because of what rsatscan requires)
# the row names in the dataframe look (and are) redundant, but that's okay, because the default of write.nbr is to ignore row names

### Write case, population, and neighbor files for SaTScan ###

setwd(dir="./SaTScan_Files/")
dir<-getwd()

write.cas(x=YawsCasesFile,filename="YawsCases",location=dir)
write.cas(x=YawsCasesFileCovars,filename="YawsCasesCovars",location=dir)
write.cas(x=YawsCasesFileCovarsConf,filename="YawsCasesCovarsConf",location=dir)

write.pop(x=AllVisitsFile,filename="AllVisits",location=dir)
write.pop(x=AllVisitsFileCovars,filename="AllVisitsCovars",location=dir)

write.nbr(x = NonEuclideanNeighbors, filename = "Lihir_Passive", location = dir)

### Run SaTScan analyses ###

# Define settings that will be used repeatedly
invisible(ss.options(reset=TRUE)) # reset satscan options to defaults
ss.options(invals=list(PrecisionCaseTimes=3, # cases are by day
                       AnalysisType=3, # retrospective space-time
                       ModelType=0, # discrete Poisson 
                       TimeAggregationUnits=3, # days are unit
                       TimeAggregationLength=7, # aggregate 7 days in a row (weekly), makes computations faster
                       OutputShapefiles="n",
                       MostLikelyClusterEachCentroidDBase="n",
                       CensusAreasReportedClustersDBase="n",
                       IncludeRelativeRisksCensusAreasDBase="n",
                       SaveSimLLRsDBase="n",
                       UseNeighborsFile="y",
                       NeighborsFilename="Lihir_Passive.nbr", # identify neighbors file
                       MaxTemporalSize=50, # 90% of time in a cluster is largest permissible size, but use 50%
                       LaunchKMLViewer="n",
                       CriteriaForReportingSecondaryClusters=0, # no geographic overlap of reported clusters
                       ReportClusterRank="y" # reports rank of clusters in Monte Carlo reps
))
ss.options(c("StartDate=2005/04/04","EndDate=2016/05/29"))


# Discrete Poisson model adjusted for age and sex (TABLE 1)

ss.options(invals=list(CaseFile="YawsCasesCovars.cas",  # identify case file
                       PopulationFile="AllVisitsCovars.pop" # identify "population" file
                       ))
write.ss.prm(location=dir,filename="YawsPoisAdj")
YawsPoisAdjAnalysis<-satscan(prmlocation=dir,prmfilename="YawsPoisAdj",verbose=T,cleanup=F,sslocation = sslocdir)


# Discrete Poisson model adjusted for age and sex, confirmed RPR+ only (TABLE S1)

ss.options(invals=list(CaseFile="YawsCasesCovarsConf.cas", # identify case file
                       PopulationFile="AllVisitsCovars.pop" # identify "population" file
))
write.ss.prm(location=dir,filename="YawsPoisAdjConf")
YawsPoisAdjConfAnalysis<-satscan(prmlocation=dir,prmfilename="YawsPoisAdjConf",verbose=T,cleanup=F,sslocation = sslocdir)


# Discrete Poisson model unadjusted for age and sex (TABLE S2)

ss.options(invals=list(CaseFile="YawsCases.cas", # identify case file
                       PopulationFile="AllVisits.pop" # identify "population" file
))
write.ss.prm(location=dir,filename="YawsPoisUnadj")
YawsPoisUnadjAnalysis<-satscan(prmlocation=dir,prmfilename="YawsPoisUnadj",verbose=T,cleanup=F,sslocation = sslocdir)


# Space-time permutation model adjusted for age and sex (TABLE S3)

ss.options(invals=list(CaseFile="YawsCasesCovars.cas", # identify case file
                       ModelType=2 # space-time permutation 
))
write.ss.prm(location=dir,filename="YawsSTPermAdj")
YawsSTPermAdjAnalysis<-satscan(prmlocation=dir,prmfilename="YawsSTPermAdj",verbose=T,cleanup=F, sslocation = sslocdir)


# Space-time permutation model unadjusted for age and sex (TABLE S4)

ss.options(invals=list(CaseFile="YawsCases.cas", # identify case file
                       ModelType=2 # space-time permutation 
))
write.ss.prm(location=dir,filename="YawsSTPermUnadj")
YawsSTPermUnadjAnalysis<-satscan(prmlocation=dir,prmfilename="YawsSTPermUnadj",verbose=T,cleanup=F,sslocation = sslocdir)

setwd("..")